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Equilibrium properties of hydrogen-helium mixtures under conditions similar to the interior of 
giant gas planets are studied by means of first principle density functional molecular dynamics 
simulations. We investigate the molecular and atomic fluid phase of hydrogen with and without the 
presence of helium for densities between p = 0.19 gem"'' and p = 0.66 gem"'' and temperatures 
from T = 500 K to T = 8000 K. Helium has a crucial influence on the ionic and electronic structure 
of the liquid. Hydrogen molecule bonds are shortened as well as strengthened which leads to more 
stable hydrogen molecules compared to pure hydrogen for the same thermodynamic conditions. The 
ab initio treatment of the mixture enables us to investigate the validity of the widely used linear 
mixing approximation. We find deviations of up to 8% in energy and volume from linear mixing at 
constant pressure in the region of molecular dissociation. 



I. INTRODUCTION 

The discovery of the first extrasolar planet in 1995i 
marked the beginning of a new era in planetary science, 
which is characterized by great improvements in observa- 
tional techniques and a rapidly expanding set of known 
extrasolar planets. Most of the about 200 known planets 
are giant gas planets in small orbits since the primary tool 
for detection, radio velocity measurement, is most sen- 
sitive for finding heavy planets that rapidly orbit their 
parent stas^. From radius measurements of transient 
extrasolar planets, we know that most of the discovered 
planets consist primarily of hydrogen and helium. There- 
fore, there is a great need for accurate equation of state 
(EOS) data for these elements under giant gas planet 
conditions^. The knowledge of equilibrium properties 
of mixtures of hydrogen and helium will help to clar- 
ify questions concerning the inner structure, origin, and 
evolution of such astrophysical objects. Open questions 
are whether or not hydrogen and helium phase-separate 
inside giant planets, if a plasma phase transitiorkSi under 
the influence of helium can be found, and if a solid rocky 
core exists in Jupiter^'^. 

The EOS of hydrogen has attracted considerable at- 
tention and a large number of models have been intro- 
duced to characterize hydrogen at high pressure and tem- 
perature. Of great use in astrophysical calculations and 
planet modeling are (free energy) models operating in 
the chemical picturaSiSiSiiSiiiiiSii^. In these models, the 
hydrogen fluid is assumed to be composed of well-defined 
chemical species like atoms, molecules, and free charged 
particles. Such methods operate in the thermodynamic 
limit and are capable of describing large parameter re- 
gions of temperature and density. Further advantages of 
the free energy models are the small computational effort 
required to calculate the EOS and the explicit knowledge 
of all the considered contributions to the EOS. Ionization 
and dissociation degrees are computed by means of mass 
action laws and are not subject to fluctuations due to 



technical issues as in simulations. Atoms and molecules 
are treated as separate elementary species instead of be- 
ing considered as bound states of electrons and nuclei. 
This implies certain approximations that limit the qual- 
ity of these approaches. At sufficiently high density, the 
definition of atoms and molecules becomes imprecise as 
the lifetime of such objects decreases rapidly with density 
and mean distances between nuclei and electrons become 
comparable to bond lengths. 

In this paper, special emphasis will be placed on testing 
the accuracy of the linear mixing approximation, which is 
often applied in free energy models to calculate the EOS 
of mixtures of different chemical species such as hydrogen 
and helium. A similar approach can be used to character- 
ize mixtures of hydrogen atoms and moleculesSiii. This 
approximation allows the calculation of thermodynamic 
variables of mixtures by a simple linear superposition of 
properties of pure substances. Linear mixing is a useful 
assumption to make if reliable experimental or theoret- 
ical data are only available for pure substances, or for 
cases where it has been shown that particles interact only 
weakly. 

To avoid the shortcomings of chemical models, first 
principle calculations can be applied. Such methods work 
in the physical picture and treat electrons and nuclei as 
elementary particles interacting via the Coulomb poten- 
tial. Quantum theory then describes the effects leading 
to the formation of atoms or molecules and their statis- 
tics. For hydrogen, there have been great efforts to study 
the equilibrium properties by means of density functional 
theory (DFT)!'"^-!*^ ", DFT-molecular dynamics (DFT- 
MD)i§^'-^S-'-?P-, DFT - hypernetted chain equation combi- 
nation (DFT-HNC)^ii^2i2^, path integral Monte Carlo 
(PIMC)24.25^ couple electron-ion Monte Carlo (QMC)i^, 
and Green's function theorjiSLSS. 

Questions addressed include the problem of the hydro- 
gen HugonioliSSiSSiSiiS^ and helium Hugoniot"^'^, the na- 
ture of the transition in hydrogen from a molecular to an 
atomic state2ii2ii2^, the melting line of hydrogen'^^, the 
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different molecular solid phasesiiai^^i^, and the atomic 
solid (metallic Wigner crystal) proposed to be found at 
very high pressures^S'^i'^. Although DFT-MD is pri- 
marily an electronic groundstate method, it can be read- 
ily applied to describe dense solid and fluid hydrogen and 
helium at conditions relevant to giant gas planets because 
the electrons in such systems are either chemically bound 
or highly degenerate. 

To our knowledge, there are only a few first princi- 
ple calculations dealing with the question of the helium 
influence on the hydrogen EOS. Klepeis et alm^ made 
predictions concerning the hydrogen-helium phase sep- 
aration but the DFT method applied in that paper is 
not suitable to treat the high temperature liquid found 
inside giant gas planets since only lattices could be dis- 
cussed. The first DFT-MD calculations, performed by 
Pfaffenzeller et al."^^, lead to more reasonable values for 
hydrogen-helium demixing. However, their simulations 
were performed using Car-Parinello MD (CP-MD). The 
high temperature region (T > 15000K) where partial ion- 
ization occurs was considered by Militzer^. 

Here we present new results concerning the hydrogen 
and hydrogen-helium EOS for conditions inside giant gas 
planets as derived from first principle DFT-MD. We use 
Born Oppenheimer MD (BO-MD) in order to ensure well 
converged electronic wavefunctions at every step. The 
density and temperature values chosen cover the region 
of molecular dissociation where we expect corrections to 
the linear mixing approximation to be most significant. 

We continue with Section HTl which contains details of 
our computational method. Results for the hydrogen 
EOS are presented in Section llH Al and compared to EOS 
data from a variety of other approaches. Furthermore, 
the ionic and electronic structures of the hydrogen fluid 
and their dependence on temperature and density are in- 
vestigated as well. The EOS and properties of hydrogen- 
helium mixtures are studied in Subsection IHI Bl The 
focus there is on understanding how the presence of he- 
lium influences the stability of hydrogen molecules and 
the electronic structure, as well as on determining excess 
mixing quantities. Finally, the validity of the linear mix- 
ing approximation is examined in Subsection IHI CI and 
Section IIVI provides a summary of our results and con- 
clusions. 



II. METHOD 

We use first principle DFT-MD within the physical pic- 
ture to describe hydrogen-helium mixtures under giant 
gas planet conditions. This means that protons as well 
as helium nuclei are treated classically. Nuclei and elec- 
trons interact via a Coulomb potential. Since T <^ Tp, 
where Tp is the Fermi-temperature, for all densities and 
temperatures found inside a typical giant gas planet, we 
employ ground state density functional theory to describe 
the electrons in the Coulomb field of the ions. The ions 
have sufficiently large mass to be treated as classical par- 



ticles and their properties described well by means of 
molecular dynamics simulations. We employ the Born- 
Oppenheimer approximation to decouple the dynamics of 
electrons and ions. The electrons thus respond instanta- 
neously to the ionic motion and the electronic wave func- 
tions are converged at every ionic time step. Compared 
to Car-Parinello MD, this reliably keeps the electrons in 
their ground state without the necessity to use an arti- 
ficial electron thermostat for systems with a small band 
gap. 

The calculations presented in this article were carried 
out with the CPMD packager's. AH MD results were 
obtained within the NVT ensemble. A Nose-Hoover 
thermostat was applied to adjust the system tempera- 
ture. The thermostat was tuned to the first vibration 
mode of the hydrogen molecule (4400 cm~^). Ah DFT- 
MD simulations were with 128 electrons in super cells 
with periodic boundary conditions and convergence tests 
were performed with larger cells. An ionic time step 
of At — 16 a. u. (la.u. = 0.0242 fs) was used through- 
out, however we found that At — 32 a. u. is already suf- 
ficient for Ts > 1.86 in the hydrogen-helium mixtures. 
{rs ~ 3/{(47rn)^/'^aB}, rg is the Wigner-Seitz radius, ?i 
the number density of electrons per unit volume). All 
simulations were run for at least 2 ps and for the calcu- 
lation of thermodynamic averages, such as pressure and 
energy, an initial timespan of at least 0.1 ps was not con- 
sidered to allow the system to equilibrate. 

The DFT calculations were performed with plane 
waves up to a cutoff energy of 35 — 50 Ha, the 
Perdew-Burke-Ernzerhof GGA approximation^ for the 
exchange-correlation energy, and F-point sampling of the 
Brillouin zone. We used local TrouUier-Martins norm 
conserving pseudopotentials^2ii&. The pseudopotentials 
were tested for transferability and for reproducing the 
bond length and groundstate energy of single hydrogen 
molecules as well as of helium dimers. 

For each density, the simulations were started at low 
temperature where the system is in a molecular phase and 
the temperature was increased in steps of AT = 500 K in 
order to avoid the premature destruction of molecules by 
temperature oscillations introduced by the thermostat. 

The electronic density of states (DOS) were calcu- 
lated for snapshots from MD simulations with the Abinit 
package^ and using a Fermi-Dirac smearing. The pre- 
sented results for DOS and bandgaps are based on mul- 
tiple snapshots for each parameter set. 

Finite size effects were tested for by carrying out sim- 
ulations with supercells containing up to 324 electrons 
(plus the required neutralizing number of protons and 
helium nuclei). For densities corresponding to 1.86 > 
Ts > 1.6 and at T = 500 K we found the changes in pres- 
sure and energy to be smaller than 2%. The convergence 
of the Brillouin zone sampling was checked by optimiz- 
ing the electronic density of several MD snapshots with 
F-point, a 2x2x2, and a 4x4x4 Monkhorst-Pack grid of 
k-points^'' for a Ne = 128 system. The Brillouin zone ap- 
pears to be sufficiently small so that deviations between 



3 



results with 1, 8, and 64 k-points are below 1% (rs = 1.6, 
T = 500 K). 

We also examined the significance of electronic exci- 
tations for the thermodynamic properties of the studied 
fluids. Snapshots from MD trajectories were taken and 
electronic states were populated according to a Fermi 
distribution corresponding to the MD temperature. For 
Ts = 2.4 and T = 7000 K the pressure was found to in- 
crease by 8%. Since the degeneracy parameter of the 
electrons increases while moving along an isentrope to 
the center of a giant gas planet, this can be considered 
an upper limit for finite temperature electronic excitation 
effects; the error for higher densities and lower tempera- 
tures will be much smaller. 
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III. RESULTS 

Here we present ab initio results for equilibrium prop- 
erties of hydrogen and hydrogen-helium mixtures in 
a density region between 0.19g-cm~'^ and 0.66g-cm~^ 
{rs = 2.4 to rs = 1.6) and for temperatures from 500 K 
to 8000 K. This parameter region includes part of the 
transition region from the molecular to the atomic fluid 
state of hydrogen and hydrogen-helium mixtures. It is, 
therefore, interesting to study not only to get insight into 
interior properties of giant gas planets but also to ex- 
amine molecular dissociation, the molecular-atomic and 
the insulator-metal transitions in hydrogen. Addition- 
ally, one can consider the influence of helium on these 
transitions and properties of mixing. 



A. Pure Hydrogen 

Figure ^ provides a summary of our hydrogen EOS 
calculations using DFT-MD. Four different pressure iso- 
chores are shown. At low density (r^ = 2.4), the pressure 
increases monotonically with temperature as the char- 
acter of the fluid changes smoothly from molecular to 
atomic. At this density, the transition is slow enough 
with temperature so that the drop in the pressure when 
molecules break (the interactions become less repulsive) 
is compensated by the increase of the kinetic contribution 
to the pressure. 

At higher density (rs < 2), the dissociation of 
molecules takes place more rapidly with increasing tem- 
perature and leads to a region of dP/dT\v < 0. At 
sufficiently high density, this effect dominates over the 
pressure increase that results from the presence of two 
atoms instead of one molecule. Furthermore, the condi- 
tion dP/dT\v < implies a negative thermal expansivity 
dV/dT\p < 0, while the fluid maintains hydrostatic sta- 
bility given by dP/dV\T > 0. 

By exhibiting a region with dP/dT\v < 0, fluid hydro- 
gen shares some properties with typical solids, where a 
new crystal structure with more efficient packing appears 



FIG. 1: Pressure-temperature relation for pure hydrogen 
along various isochores in DFT-MD (this work) and accord- 
ing to Saumon and Chabrier (SC)i-. The isochore of SC for 
rs = 1.6 lies out of the range of the P axis. Errors for the 
DFT-MD simulations are of the order of the symbols. 
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FIG. 2: Pair correlation functions g{r) for pure hydrogen 
at different temperatures (color coded) and different den- 
sities: solid, dash-dotted, and dashed lines, respectively: 
rs = 1.6 (0.66 g ■ cm"^), = 1.86 (0.42 g • cm"^), = 
2.4(0.20g-cm-3). 



and the pressure is lowered at fixed volume. In solid hy- 
drogen different transition pressures to an atomic solid 
have been predicted, above SOOGPa^i*^ and =^ 1.3. 
This is consistent with the observed shift of the region 
of dP/dT\v < to lower temperatures, as the density 
is increased. It indicates a density effect on dissocia- 
tion as less and less thermal energy is needed to break 
up the molecular bonds. At even higher densities than 
shown here, the bond length is equal or less than the 
mean particle distance. In this regime the interaction of 
molecules and atoms becomes too strong and pressure 
dissociation/ionization occurs. 

Figure[21shows the dependence of pair correlation func- 
tions (pure hydrogen at r^. = 2.4) on the temperature. 
The first peak at r w 1.40 a.u. indicates the existence 
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of hydrogen molecules With increasing temperature, the 
height of the first peak is reduced as molecules dissoci- 
ate. An analogue behavior can be observed by plotting 
the changes in g{r) with density. A strong decrease of 
the first peak with increasing density and thus a signif- 
icant lower fraction of molecules at higher densities is 
revealed. In addition to the less pronounced first maxi- 
mum, an overall weakening of the short range order can 
be observed with increasing temperature. 

A more quantifiable picture of the described effects can 
be obtained by plotting the degree of dissociation. 
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as a function of temperature (Fig. ISJ. Here is the 
average number of hydrogen molecules at the given den- 
sity and temperature conditions. Nh the total number of 
hydrogen nuclei irrespectively of the dissociation state. 

At higher density, fewer molecules are present at the 
same temperature as a result of pressure dissociation. 
While at low density, the dissociation proceeds gradually 
with temperature, the curves for rg < 1.75 show a rapid 
drop around 2500 K, which is related to the dP/dT\v < 
region. 

The dissociation degree and the binary distribution 
function are nevertheless not sufficient to draw a com- 
plete picture of the structure and dynamics in fluid hy- 
drogen. The lifetime of the molecules must also be taken 
into account. Figure |3| shows that at rg = 1.75, for exam- 
ple, even though on average more than 50% of the protons 
are found in paired states, the lifetime of these pairs is 
short (less than 2 H2-vibrations on average); there is a 
continuous formation and destruction of pairs of hydro- 
gen atoms. It is therefore imprecise to classify the fluid 
as either molecular-atomic or pure atomic as there is no 
unique criterion for a molecule. However, the results for 
the EOS obtained by our simulations do not depend on 
the number of molecules or atoms but only on tempera- 
ture and density. 

In addition to the dissociation of hydrogen molecules, 
of interest are the changes in the electronic structure 
taking place in the same parameter region. As the sys- 
tem becomes denser or the temperature is raised (still 
T <C Tp), the interactions between the molecules in 
the fluid become stronger, and the formerly well-bound 
electrons become delocalized. This is associated with a 
strong increase in the electrical conductivity and is usu- 
ally referred to as metallization^. The effect can be 
seen in the electronic density of states (DOS), namely 
the band gap, as shown in Fig. ^ We calculated the 
Kohn-Sham eigenvalues in GGA for several snapshots 
and estimated the bandgap in fluid hydrogen along the 
MD trajectory. For — 1.6, a closing of the bandgap 
can be observed around T = 2000iir. This means that a 
metallic-like state may have been formed. Furthermore, 
as indicated by the red lines in Fig. 21 the degree of dis- 
sociation incorporating lifetime effects decreases strongly 
around the same temperature. To distinguish this degree 
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FIG. 3: Dissociation degree of the hydrogen molecules in pure 
hydrogen. The dissociation degree obtained by simply count- 
ing all pairs of hydrogen atoms with distance shorter than 
rcut ~ 1.8 a.u. is plotted with dotted lines. Full lines take 
into account the lifetime of these pairs as well (10 H2 vibra- 
tions at least to be counted as molecule). 
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FIG. 4: Bandgap calculated in GGA for fluid hydrogen and 
hydrogen-helium mixture (blue curves) as well as dissociation 
degree for the two systems (red curves). The density corre- 
sponds to ra — 1.6. The dissociation degree is determined by 
taking into account lifetime effects. 



of dissociation we do not consider proton pairs with less 
than 10 H2 vibrations (t = 10 x 7.6 fs) as being molecules. 
The closing of the bandgap and dissociation of hydrogen 
molecules happen at the same time. 

It is well known that the GGA underestimates the 
bandgap. More sophisticated calculationsi^ (for a H2 
sohd) give a bandgap of ~ 0.6 Ha at T = 300 K. The 
value we obtain is 4 times smaller. Hood and Galli^ 
compared values obtained with DFT (GGA) to quan- 
tum Monte Carlo (QMC) gaps for liquid deuterium. For 
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T = 3000 K and = 1.6, they found the QMC gap 
twice as large as the DFT gap. The actual temperature 
of metallization is thus somewhat higher. To improve 
the description of the electronic properties of the fluid, 
one needs to use a more accurate method than DFT- 
GGA, which is beyond the scope of this article. Still, it 
is worth noting that we find a continuous transition from 
an insulating to a conducting state, as determined by the 
closing of the gap. We do not observe molecules in the 
conducting phase as found by Weir et alP^ or Johnson et 

While there is a general agreement about dissociation, 
ionic, and electronic structural changes throughout var- 
ious papers, this agreement is only qualitative. As can 
be seen in Fig. ^ different methods give very different 
results for the EOS of dense fluid hydrogen. Whereas for 
the lowest density shown in Fig. ^ the agreement be- 
tween the free energy model of Saumon and Chabrier^s^ 
and our results is reasonable, deviations up to 20% (at 
Ts = 1.86) and even 24% (at Ts — 1.75 and above) can 
be found for higher densities. The free energy model 
overestimates the pressure considerably. The degrees of 
dissociation calculated with the free energy method show 
a significantly higher fraction of molecules than we find 
in our simulations. However, even a linear extrapolation 
of our pressure results of the molecular phase to higher 
temperatures (a linear scaling very similar to the one at 
Ts = 2.4 is assumed) only reduces the discrepancy with 
the result of Saumon and Chabrier but can not elimi- 
nate the difference completely. Deviations of this order 
may signify completely different physics inside giant gas 
planets and it is of great importance to discuss the dis- 
crepancies. 

Comparisons with different first principle calculations 
can help resolve this issue, as agreement between different 
independent ab initio methods would be a strong indica- 
tion for correct results. Figure |S1 provides such a com- 
parison. In addition to our results, PIMC data™, wave 
packet molecular dynamic (WPMD) rcsultsSLSMi and 
older DFT-MD points^ are shown. Furthermore, the 
isochores of two different models in the chemical picture 
are added: fluid variational theory (FVT)5£ as well as the 
linear mixing (LM) modeU^. They start from a mixture 
of atoms and molecules and their (Lennard-Jones type) 
interactions to minimize the free energy with respect to 
the fraction of the constituents. The isochore provided 
by WPMD deviates from the other ones by more than a 
factor of two at lower temperatures and by 25% at the 
highest temperatures shown here. PIMC is not a ground 
state method and is more capable of determining the EOS 
at higher temperatures. Information about fluid hydro- 
gen or even about dissociation of hydrogen molecules can 
not be obtained. For higher temperatures, PIMC results 
lie between DFT-MD and FVT data. Better agreement 
is achieved between the two DFT-MD methods, FVT and 
LM. In the region with temperatures less than 10000 K 
some features of the curves are still different. The first 
principle simulations show a region with a reduced or 
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FIG. 5: Pressure-temperature relation for a single isochore 
of hydrogen with different methods for = 2. Our DFT- 
MD results are shown in blue with circles. FVT by Juranek 
et alm^ (sohd black hne), PIMC by Militzer et ali^ (dot- 
ted line, triangles), WPMD by Knaup et ali~ (dashed line, 
squares), DFT-MD by Collins et al^ (dashed, diamonds), 
LM by RosE^" (long-short dashed), and this works DFT-MD 
results (blue, dots). The black dot indicates the pressure of a 
Hz-solid at r = 300 K^S. The error bars of simulation results 
are within the size of the markers. 



even slightly negative slope (smooth transition from a 
purely molecular to an atomic fluid) around 5000 K. Such 
a behavior is absent in the FVT model. The LM method, 
on the other hand, shows a similar feature but on a much 
wider temperature range. Thus, there is no unique pic- 
ture of the EOS of dense fluid hydrogen. Deviations be- 
tween the results of different first principle methods are 
related to the treatment of the electrons. Inconsistencies 
between the chemical picture, as a basis for FVT and LM, 
and the physical picture, as a foundation of DFT-MD or 
PIMC, contribute to the non-unique description. How- 
ever that may be, DFT-MD provides reproducible results 
(the two DFT-MD studies were performed with different 
codes) and FVT seems to be reliable in the molecular 
fluid phase. 

Lastly, we consider more closely the intermediate den- 
sity range at a temperature typical for the interior of 
Jupiter. Figure El shows such an isotherm for T = 
5000 K. Our result gives the lowest pressure. Saumon 
and Chabrier's EOS predicts a pressure up to 20% higher 
as stated above. The blue curve, as well as the green 
one, show a plasma phase transition (PPT). According to 
Saumon and Chabrier the PPT is expected at a density of 
approximately p = 1 g-cm^^ (p ^ 200 GPa). The FVT 
results withJi or without^*' extension to ionized plasmas 
gives an isotherm right between the two above mentioned 
results. These methods predict a PPT at slightly smaller 
density of p = O.Sg-cm-^ and P = 100 GPa). At both 
densities, neither we nor Weir et alm^ have evidence for 
such a PPT. Instead, a continuous transition from a 
molecular to an atomic state takes place at lower temper- 
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FIG. 6: The 5000 K pressure isotherm of hydrogen com- 
puted with various methods: this works DFT-MD (black 
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Saumon and Chabrier (blue solid)". HNC-MAL-- gives simi- 
lar results as FVT. 



atures. It is nevertheless remarkable that the inclusion 
of ionization reduces the pressure and gives better agree- 
ment. From our simulation we have no information about 
possible intermediate ionized states of hydrogen. 



B. Hydrogen-Helium Mixtures 

So far we have studied pure hydrogen under giant gas 
planet conditions. A further degree of freedom is added 
if one considers a mixture of hydrogen and helium. He- 
lium, even in small fractions, changes the EOS signifi- 
cantly. Helium has an influence on the formation and 
dissociation of hydrogen molecules, and it changes the 
ionic structure of the liquid as well as the electronic prop- 
erties. The transition from a molecular state into an 
atomic state may be displaced or its character changed. 
Further, hydrogen and helium have been predicted to 
phase-separate in giant planet interiors. ^iSSiS. 

Models in the chemical picture use the linear mixing 
rule to add hydrogen and helium portions to the EOS®^. 
Contributions from the entropy of mixing are ignored 
and all the interactions between the two subsystems are 
left out. First principle calculations include all these ef- 
fects since a mixture of the two fluids can be simulated 
directly. The demixing line was calculated by classical 
Monte Carlo^, by ground state DFT calculations^^, and 
by Car Parinello MD^^ 

Here, we primarily present results for a hydrogen- 
helium mixture at a mixing ratio oix = 0.5. The mixing 
ratio is defined as 
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where Njj and Njje are the number of hydrogen and he- 
lium nuclei per unit volume. This definition weights the 
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FIG. 7: Pressure isotherms of a mixture of hydrogen and he- 
lium (x=0.5) for different temperatures. High density limiting 
results for hydrogen and an electron gas are shown addition- 
ally. 



species according to the number of electrons that they 
contribute to the system. The corresponding Wigner- 
Seitz radius is computed from the total number of elec- 
trons. For many simulations, a mixing ratio of a; = 0.5 
was chosen so that large interaction effects between the 
two species could be observed. 

Figure [7| shows the pressure for a number of isotherms 
for a hydrogen-helium mixture. The maximum density 
shown here corresponds approximately to conditions in 
the center of Jupiter (r^ = 0.9, p = 3.6gcm~'^). It 
is demonstrated that temperature is not important for 
higher densities, since all the isotherms merge into the 
one with the lowest temperature. At the highest densities 
shown here the temperature contribution of the ions to 
the pressure is approximately 5%. The ions are strongly 
coupled and their interaction contribution to the EOS is 
of the order of 30%. The rest of the deviation from the 
ideal degenerate system is given by nonidealities in the 
electron gas and interactions between electrons and ions. 
For even higher densities the electronic contributions will 
become even more important since they rise with density 
as n^/^ and thus faster than any other contribution. 

The role of helium for the EOS of the mixture can 
be studied in Fig. The pressure is slightly lowered 
over the whole temperature range, but more important 
is the fact that the region with negative dP/dT\v has 
vanished (for x = 0.5). The additional curve for Jupiter's 
hydrogen-helium mixing ratio (x=0.14) shows an inter- 
mediate step where the region with negative slope of the 
pressure does still exist but is reduced in markedness and 
in the temperature range of existence. More helium in the 
mixture means a shift of the negative slope to higher tem- 
peratures and a decrease of the depth of the minimum. 
For very high temperatures, the EOS of pure hydrogen 
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FIG. 8: Pressure isochores for a mixture of hydrogen and he- 
lium (x=0.5, solid lines) and for pure hydrogen (dotted lines) 
at two different densities. For = 1.6 an additional curve 
for Jupiter's helium ratio (x=0.14) was added (back dashed 
line) . 




I 1 I 1 I 1 I 1 I 1 I 1 L 

1000 2000 3000 4000 5000 6000 7000 



T[K] 

FIG. 9: Comparison of dissociation degrees (incorporating 
life time effects) in pure hydrogen and in hydrogen- helium 
(x=0.5). The solid lines are nonlinear fits to the data points. 



and mixtures show a very similar behavior. Hehum has 
a more significant influence on the pressure at lower tem- 
peratures. 

These two different regions can clearly be assigned to 
different dissociation regimes (see Fig. I^l. The transi- 
tion from a molecular phase to an atomic phase, while 
still smooth, takes place at lower temperature and over 
a shorter range of temperature in hydrogen than in the 
mixture. This relatively rapid change in the microstruc- 
ture of the fluid is the reason for the drop in the pressure. 



The vanishing of the molecules in hydrogen and their ex- 
tended existence in the mixture can be confirmed with 
the help of pair correlation functions in Fig. ^| The 
molecular peak (first peak) drops considerably faster in 
pure hydrogen. This behavior and the higher peak in the 
mixture give clear evidence for molecules at the high- 
est temperatures shown here. If we take the position 
of the first maximum as a measure for the mean bond 
length of the hydrogen molecule, we obtain a value of 
(d) = 1.37 ao for pure hydrogen. In the mixture (with 
ratio of a; = 0.5) this value changes to (d) — 1.29 ao which 
means a shortening of the bond by 6%. The same can 
be obtained by means of nearest neighbor distributions 
as in Fig. ^] The first neighbor distribution considers 
the nearest neighbor only and effects of particles farther 
away are removed from the curve. Bond lengths obtained 
from Fig. ^Ja-re slightly bigger. In addition, a shift of 
the bond length in hydrogen from 1000 K to 2000 K is 
revealed. The reduction of the bond length of 6% is con- 
firmed. The latter value is in rather good agreement 
with data by Pfaffenzeller et alm^. The same conclusion 
is derived when comparing pair correlation functions for 
pure hydrogen and hydrogen-helium mixtures at constant 
pressure instead at constant electronic density. 

Figure |51 describes the role of the (electronic) den- 
sity during the process of dissociation. Helium stabi- 
lizes the molecules. This is due to the higher charge of 
Z = 2 of the helium nuclei. The intramolecular bonds de- 
pend strongly on the electronic behavior and the space 
available. If electronic wavefunctions between different 
molecules start to overlap, in other words, when Fermi 
statistic becomes important for the electrons of the sys- 
tem as a whole, and when the distances between the 
particles become so small that interactions between the 
molecules are no longer weak, the electrons are forced to 
delocalize to obey the Pauli exclusion principle and bond- 
ing becomes impossible. Helium under giant gas planet 
conditions, in atomic form, binds two electrons closely. 
The rest of the hydrogen atoms and electrons are af- 
fected less by density and temperature and the molecules 
remain stable over a wider range. The helium influence 
is thus two-part. First, its stronger Coulomb attraction 
binds electrons. Second, (as a consequence) it influences 
the many particle state of the electrons. 

A phase diagram for molecular and atomic hydrogen 
and hydrogen-helium is provided in Fig. 1121 Here, pa- 
rameter regions for purely molecular, purely atomic as 
well as intermediate phases are shown. The diagram 
shows the increasing differences between hydrogen and 
the mixture with increasing pressure (density) and the 
huge differences (especially at high pressure) in the rate 
of the transition from a molecular to an atomic state. 
The transition region (from 95% to 5%) has nearly the 
same size for small pressures in pure hydrogen and in the 
mixture. At the other end of the pressure scale hydrogen 
changes from molecular to atomic over only 1500 K. In 
the mixture the changes are more moderate. At the pres- 
sures considered here, pressure dissociation is suppressed 
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FIG. 10: Pair correlation lunctions for hydrogen (dashed 
lines) and hydrogen-helium (solid lines) at different tempera- 
tures. The order of the curves from top to bottom is the same 
as in the legend. The vertical thin dashed lines indicate the 
location of the first peak for hydrogen and hydrogen-helium, 
respectively. 
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since the slope of the lines with constant dissociation de- 
gree is small for higher pressures. A good zeroth approxi- 
mation for the behavior of the mixture is given by taking 
into account the hydrogen density only for the dissoci- 
ation process. In this way, = 1.6 for a 50% mixture 
of hydrogen and helium would correspond to Vg = 2.02 
in pure hydrogen. This estimate works quite well for the 
5% line, for instance. 

Band structure and electronic density of states are af- 
fected also by helium. The inner regions of Jupiter are 
believed to be made of helium rich metallic hydrogen-*. 
The conditions under which the mixture becomes metal- 
lic strongly depends on the amount of helium. A com- 
parison of the (GGA) bandgaps in pure hydrogen and a 
hydrogen-helium mixture (x = 0.5) is shown in Fig. 0] 



8000 



7000 



6000 



5000 



4000 



3000 



2000 




20 30 



40 50 60 
P [GPa] 



80 90 100 



FIG. 12: Temperature-pressure-plane showing estimated lines 
of constant dissociation degree in pure hydrogen and in a 
hydrogen-helium mixture (x=0.5). The percentages give the 
fraction of hydrogen atoms bound in molecules. 
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FIG. 13: Electronic density of states for pure fluid hydrogen, 
pure fluid hehum, and for a mixture of hydrogen and helium 
(x=0.5) at fixed T = 500K and = 1.86. The pressure is 
approximately 30 GPa. Curves shifted to agree at the Fermi 
energy (0 Ha) which is taken in the middle between HOMO 
and LUMO. 



And while the bandgap in pure hydrogen goes to zero 
at relatively low temperatures, the gap in the mixture 
remains open over the whole temperature region shown. 
This can be traced back to the charge of the helium nu- 
cleus which shifts part of the Kohn-Sham eigenvalues to 
lower energies and thus increases the gap. The change 
in the electronic DOS from a helium system to a mix- 
ture to a pure hydrogen system is shown in Fig. 1131 The 
black peak indicates atoms in fluid helium. The red curve 
shows mainly molecules in hydrogen and there are peaks 
resulting from intermolecular interaction, too. The blue 
curve for the mixture is a superposition of the ones for 
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FIG. 14: Electronic density of states for pure fluid hydrogen 
and for a mixture of hydrogen and helium (x=0.5) as function 
of temperature. The pressure is approximately 100 GPa. 

pure systems, namely a helium peak to the left and a hy- 
drogen molecule peak on the right hand side. The bands 
on the right hand side are empty since the curves are 
normalized so that the Fermi energy is at Ha. The po- 
sition of the edges of the bands at positive and negative 
energy strongly depends on the amount of helium in the 
fluid and therefore the width of the bandgap depends on 
the helium amount. 

The dependence of the DOS on temperature is demon- 
strated in Fig. El We observe a similar effect as Scan- 
dolo dicl2^ although the transition from a mainly molecu- 
lar liquid to an atomic fluid is smooth in our calculations. 
Due to the larger initial bandgap in the hydrogen-helium 
mixture and due to the effect of the helium described 
above, the bandgap for this systems remains, even at the 
highest temperature shown. Conversely, the gap in the 
pure system has completely closed at T = 4000 K. 

C. Thermodynamic Properties of Mixtures 

We will test the validity of a common approximation 
used in determining EOS's of mixtures. Most of the effort 
has been put into the determination of the EOS of pure 
systems. By ignoring the exact nature of the interaction 
between the pure phases, one can construct the EOS of 
any mixture of the original phases with the help of the 
linear mixing (LM) approximation 

n.M = (i-^)ni+^nie' (3) 

with X according to Eq. (0) being the fraction of helium 
in the mixture and Y being a thermodynamic variable 
such as volume, pressure, or internal energy. For the free 
energy, an additional term describing the entropy of mix- 
ing must be includedi^. Linear mixing may be performed 
at constant chemical potential, at constant volume, or at 
constant pressure. For calculations of the internal struc- 
ture of giant gas planets, mixing under constant pressure 



is the most important. In all cases, the deviation from 
LM can be calculated as 

A>niix(^) = ^(-^)->LM(^)' (4) 

where Y{x) is the value obtained by DFT-MD for mix- 
ing fraction x and Y^j^ is the LM value computed from 
independent simulations results for pure hydrogen and 
pure helium. 

In this way it is assumed that the potential between 
particles of two different species can be written as an 
arithmetic average over the interactions in the pure 
systemgi^. This of course works for weak correlations 
only. The advantage is that one does not need to know 
exactly the interaction between, e.g., hydrogen and he- 
lium. This would be necessary for models in the chemical 
picture. Therefore, this approximation is used mainly by 
chemical models for the description of hydrogen-helium 
mixturesiSiii^ or even mixtures of atoms and molecules 
of hydrogeni^. The error introduced is difficult to quan- 

First principle calculations are able to verify the as- 
sumptions for LM and the validity of the approxima- 
tion since these methods rely on the more fundamental 
Coulomb interaction, do not need to assume different in- 
teraction potentials between different species, and can 
thus simulate mixtures directly. There have been some 
investigations concerning the validity of LM by classical 
MC and integral equation techniques for classical binary 
liquids not including moleculeaSSiSLSSiSS. In these cases, 
the deviations from LM found are of the order of 1% and 
below. 

The first PIMC calculations for hydrogen-helium mix- 
tures found deviations from LM at constant volume of 
up to 12% for temperatures between 15000 and 60000 K 
and giant gas planet densities {rg = LSG)"*"*. This gives 
reason to expect that linear mixing might give a slightly 
falsified picture of the EOS of hydrogen-helium at lower 
temperatures, too. Again, since the conditions for phase 
separation depend strongly on small changes in the EOS 
(and thus from deviations from LM) it is crucial to inves- 
tigate LM£&. The advantage of first principle calculations 
is the correct treatment of the degenerate electrons and 
bound states, which is missing in the classical simula- 
tions. 

The error introduced due to LM at constant volume 
(constant electronic density) as observed within DFT- 
MD is plotted in Fig. E] Similar to the findings of 
other authorsMiSS, the error is positive. As expected, the 
deviation in the pressure from the LM value is biggest 
for X = 0.5. Furthermore, increasing temperature causes 
an increase in the LM error of up to 12% for the highest 
temperature plotted in Fig. 1151 The deviation from LM 
for a Jupiter like mixing ratio of a: « 0.14 ranges from 
around zero (T 500 K) up to 10%. The temperature 
inside Jupiter for this density is according to Saumon and 
Chabrieria of the order of 5000 K. This means that one 
can expect a deviation of approximately 10% of the true 
EOS from the one calculated with LM. 
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FIG. 15: Mixing error in the pressure due to the hnear mixing 
approximation at constant volume for various temperatures 
as function of the mixing ratio. The electronic density is 
ra = 1.86. This corresponds to q = 0.42g-cm~'^ for pure H 
and Q = 1.66g-cm~'^ for pure He (pressure between 10 GPa 
and 40 GPa). The symbols represent calculated values, the 
lines were obtained with a polynomial fit of fourth order. 
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FIG. 16: Mixing error in the pressure due to the linear mix- 
ing approximation at constant volume for various densities 
as function of temperature. The mixing ratio is a; = 0.5. 
The rs — 2.4 data points suffer from simulation noise and 
the curve was obtained by least square fitting a third degree 
polynomial. 

The dependence of the error maximum at x = 0.5 from 
density and temperature is shown in Fig. 1161 The curve 
for the smallest density shown here (r^ = 2.4) gives rea- 
son to conclude that LM is a good approximation for the 
pure molecular phase of hydrogen and hehum as found 
at this density over a wide temperature range. With in- 



FIG. 17: Mixing error in the volume due to the linear mixing 
approximation at constant pressure for hydrogen-helium for 
various densities and two different mixing ratios of x = 0.5 
(solid lines) and a; = 0.14 (dashed line) as function of temper- 
ature. 



creasing density the deviations from LM start to grow as 
well and corrections to the pressure become significant for 
smaller temperatures. 5% error is reached around 3000 K 
for Ts = 1.86, around 2500 K for rg = 1.75, and at approx- 
imately 1250 K for Ts — 1.6. The maximum of APmix/P 
is located at a slightly higher temperature than the tran- 
sition from a pure molecular to a mainly atomic phase 
in pure hydrogen. The linear mixing rule transfers the 
behavior of pure hydrogen in an incorrect way into the 
mixture. This causes these deviations of up to approxi- 
mately 15%. In addition, it is shown that linear mixing is 
not a good approximation for hydrogen-helium systems 
containing atoms and molecules. For higher tempera- 
tures the deviation from linear mixing declines although 
in the considered range it does not reach values below 
5% again. 

A similar statement can be made for mixing at con- 
stant pressure as shown in Fig llTI The same features as 
in Fig. 1161 can be observed. The maximum of the mixing 
error is shifted to lower temperatures for higher densities. 
However, the error in the volume introduced by LM is 
slightly smaller than the one in the pressure. Comparing 
curves with different mixing ratios at constant density in 
Figs. 1171 and il8l it is obvious that the maximum in the 
error is reached at lower temperatures for smaller mixing 
ratios x. This is in agreement with the temperature shift 
of molecular dissociation as function of the helium ratio 
in the system. Less obvious is the actual absolute value 
for the deviation from linear mixing. Whereas the error 
in the volume never exceeds 5%, the energy is much more 
sensitive to deviations from LM with an error of up to 
9%. 

All, figures [TBlll7l andll8lshow that LM can be consid- 
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FIG. 18: Mixing error in the energy due to the linear mixing 
approximation at constant pressure as function of the temper- 
ature for hydrogen-helium mixtures at a density of Ts — 1.86. 
The mixing ratios are x = 0.5 and x = 0.14 (Jupiter). 

ered a good approximation only for systems at low tem- 
peratures (T < lOOOX) and at very high temperatures 
T > SOOOif for the densities presented here. In the case 
of liquid hydrogen-helium mixtures, these systems con- 
sist of weakly interacting hydrogen molecules and helium 
atoms (low T) or of weakly interacting hydrogen and he- 
lium atoms (high T) . More complicated situations where 
atoms and molecules of hydrogen and helium are involved 
and interactions between them are non negligible require 
a better description than LM since helium has a signifi- 
cant influence on the dissociation degree and the mixture 
cannot be considered to be a composition of two fluids. 

The figures presented here could suggest that LM is 
a rather good approximation for all of the higher tem- 
peratures. This is by no means the case. As stated 
before, LM works well only for nearly ideal systems. 
When the temperature becomes too high so that the 
gas of hydrogen and helium atoms experiences ioniza- 
tion (this can be accomplished by increasing the density 
as well) and a partially ionized plasma is created, the 
merely short ranged interatomic interactions are replaced 
by long range Coulomb forces and nonideality contribu- 
tions to the EOS become very important again. In this 
regime, LM breaks down again, as was demonstrated by 
Militzer by means of PIMG^. 

IV. SUMMARY 

We use first principle DFT-MD simulations to study 
equilibrium properties of hydrogen and hydrogen- helium 



mixtures under extreme conditions. The results obtained 
are relevant for the modeling of giant gas planets and for 
the principle understanding of the EOS of fluid hydrogen 
and hydrogen- helium mixtures. 

Our results for pure hydrogen show a smooth transition 
from a molecular to an atomic state which is accompanied 
by a transition from an insulating to a metallic like state. 
In the transition region, we find a negative temperature 
derivative of the pressure. The results for the hydrogen 
EOS show deviations from widely used chemical models 
(up to 20%). In particular, the point of dissociation for 
the molecules is obtained at much lower temperatures 
than in chemical models. We find satisfying agreement 
with previous DFT-MD simulations only. 

In particular, we demonstrate the influence of helium 
on hydrogen molecules. The presence of helium results 
into more stable molecules and an altered transition from 
a molecular to an atomic fluid state. Helium reduces the 
negative slope of the pressure isochores in the transition 
region. The bond length of the hydrogen molecules is 
shortened by 6% for x = 0.5. As a result, the degree 
of dissociation is lowered and the electronic bandgap is 
increased The effect of helium is found to be more impor- 
tant for higher densities where the stronger localization 
of the electrons prevents degeneracy effects for the elec- 
trons to become dominant. 

Our analysis of the mixing properties for a a; = 0.5 
mixture of hydrogen and helium shows that the correc- 
tions to the linear mixing approximation are significant. 
Maximum EOS corrections of 15% were found for mixing 
at constant volume and 8% for mixing at constant pres- 
sure. For Jupiter like conditions, corrections up to 5% 
were obtained. 

The presented results and forthcoming work should 
help to clarify long standing questions concerning the 
formation process of giant gas planets, help restrict the 
core size of Jupiter, and allow one to make predictions 
for the hydrogen-helium phase separation. 
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